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In this paper a general discussion of the resonance cases in an axially-symmetric potential field is pre- 
sented, when the unperturbed frequencies in the radial and z direction have a rational ratio. The general 
form of the third integral is not valid in these cases because of the appearance of divisors of the form 
(m 2 P— n*Q), which become zero in the resonance cases. However, a new isolating integral of the unperturbed 
case is available, and this can be used to construct a third integral in the form of a power series and eliminate 
all secular terms. Three cases are distinguished, (a) m+n> 4, (0) m-fw = 4, and (y) ?n+n<4. In the first 
case the orbits are rather similar to those of the general irrational case. In the third case the orbits show a 
quite peculiar character, which, however, can be explained rather accurately by a first-order theory of the 
third integral. Numerical integrations were made for the cases F = 16Q, 4F=9Q, and F = 4Q. The third 
integral, given in first- or second-order approximation, is rather well conserved. Case /3 and the cases of 
small divisors, when is near zero but not equal to zero, are discussed in Paper II. 


1. GENERAL THEORY 

I N a previous paper (Contopoulos, 1960, referred to 
hereafter as paper A) a formal third integral of 
motion is given in the axially symmetric field 

1 / C2 \ 

U — — ( +P?+Qz>-2b& , (1) 

2 \ r 2 / 

where r, 0, z are cylindrical coordinates, £=r— r 0 (r 0 , 
zc=0) is the initial point on the plane of symmetry, C is 
the angular momentum along the z axis, and P, Q , and 
b are constants. 

Then if R, Z are the projections of the velocity v on 
the axes r and z, the energy integral is 

F= i(P?+R 2 +Qz 2 +Z 2 - 2 b&) = W- (2) 


then the frequencies of the oscillations along the axes r 
and z in the unperturbed case (£ = 0) are commensurable 
and the orbits are closed. We shall see that in these 
cases also there is a formal third integral of motion. 

If $ is an integral of motion then the Poisson bracket 
(<i>,F) is zero. If we write F = Fo+bFi and <!>=<I>o+M>i 
-\-b 2 $> 2 ~\ we have 

(*o,Po) = 0, 

(*ofd+(*ifFo) = 0, 

and generally 

(^,P 1 )+(^ + i,Po) = 0, 

or 


d$?+i 

R +Z P£ Qz = - 

d( dz dR dZ 


(4) 


The third integral is found in the form of a formal series 

b 

4>=f(P£ 2 +F 2 )+ — [(P — 2Q) £z 2 — 2 £Z 2 + 2FzZ] 

4 Q-P 


— — r_ 

'-p)(o-p)L 


(Q— P)z* (2P-\-Q)R-Z 2 


(4Q-P)(Q-P)i 2 2 PQ 

(: 2P-5Q)?Z 2 (4 P-Q)RW 


2Q 


2P 


+IQI 2 z 2 +6|X2ZjH 

Ro 2 b 2 (2P+Q)R<?Zo 2 


2 0 IQ-P)(Q-P)2PQ 


(3) 


It can be seen that each term has a number of divisors 
of the form (m^P—n^Q). There is an infinite number of 
such divisors, so that if P^/Q^=n/m= rational, one such 
divisor becomes zero, and the above series is no longer 
valid. We call these cases “resonance cases” because 


Thus 4\+i is found by solving (4) if 4> 0 , 4>r • -4>, are 
known already. The corresponding system to Eq. (4) is 


dz dR dZ d $ v + 1 


= dT, (5) 


R Z -Pi -Qz — 

where T is an auxiliary variable. 

The solution of this system is 

(2*„)» (2V 0 )» 

£= smPi(T-To), 2 = sin Q'T, 


Pi 


Q k 


R= (24> 0 )i cos P*(T- To), Z= (2 V 0 )» cos Q»r, 


( 6 ) 


where 4>o, F 0 and T 0 are 3 integrals of the system (5), 
and i>o-\-Vo=Fo. Then 


— — / (*,fi)dT+f(*o,Vo,T 0 ), 


(7) 


where / is an arbitrary function of 4> 0; V 0 and P 0 - 
In order to do the integration (7) £, R , z, Z are re- 
placed by their values (6) and after the integration the 
sines and cosines are given again as functions of £, F, 
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z, Z. This is possible whenever (4>„,Pi) does not contain 
a constant term, but only trigonometric terms of the 
form 

sin 

[ffiPKr-r 0 )-nQ*r]. (8) 

cos 

If, however, a constant term appears, then the inte- 
gration gives a secular term in T. It has been proved 
in Paper A that secular terms do not appear if P*/0* 
= irrational; but if P^/Q^ = n/m, then there appear 
terms of the form 

sin sin 

[mPKT- To) - w0*P] = ( - mP*T 0 ) = const. 

cos cos 

in some (4>„, F i), that give rise to secular terms in 

In order to avoid such terms, Whittaker (1916, 1937) 
has used two methods, that in our notations can be 
described as follows: 

(a) We may choose the function /(4>o,Fo,Po) at each 
step in such a way that no secular terms appear in 
higher-order terms. 

(b) The above method does not work if the critical 
terms appear in 4>i, i.e., if 40— P=0. In this case we 
assume first that P^40 and multiply both members 
of Eq. (3) by (40— P); in this way we have the integral 

(40— P)4>= (40— P)4>o 

+^[(P— 20)£z 2 — 2£Z 2 +2P0Z]+* • 

If now 40— P becomes zero, we get 

$'=(P-20)?2 2 -2?Z 2 +22PZ+* • •. (9) 

Whittaker gave two examples to illustrate his 
methods, but he did not prove that these methods 
always work, i.e., that we can always avoid the secular 
terms. 

An indirect proof that in the resonance cases also 
there are formal integrals of motion was given by Cherry 
(1914). Cherry also gave some examples where new 
integrals are found. 

We shall apply presently a general method, that may 
be considered as a combination of the two methods of 
Whittaker, to prove that in the case of the potential 
field (1) we always have a formal third integral of 
motion, whenever P*/Q*=n/m. 

If P^/ty^n/m we have three isolating integrals of 
the system 

d£ dz dR dZ 

-=-= = . ( 10 ) 

R Z -P£ -Qz 

namely $0? V q and 

5 0 = (2<J> 0 ) im (2Fo) in smimP'To) 

= (23>o) im (2Fo) in sinlnQ'T-mPHT-To)']. 


In fact using the solutions (6) we get 
sin£wP 4 (T— T o)] 

1 

[C m iP*ZR m - 1 -C m i P*?R m - 3 -\ ], 

(2$ 0 ) !m 

1 

cos[wP*(P-Po)]= lR m -C m 2PeR m ~ 2 +-- •]. 

(2<I>o) im 

1 

sin[«0 4 P]= lC n 1 QhZ n - 1 -C n zQh*Z n ~ 3 +-- •], 

(2F 0 ) 4n 

1 

cos [«g 4 P]= [Z n — C» 2 Qz 2 Z"- 2 + • • •], 

(2F 0 ) 4n 

hence 

lim ] 

s 0 = E E (-i )«+« 2 

< 71=0 < 72*=0 

XC m 2n C\q 2+ i(P^) in R m - 2qi (Q i z) Mli - l Z n - 2t r 1 

H(>»-1)] tin] 

- E Z(-i ) SI+82 

9 1=0 92=0 

XC-2 91 + iC w 2 52 (P^) 2 ^ +1 P- 2 ^- 1 (0^) 2 ^Z«- 2 ^ (11) 

where the symbol [A~] means the biggest integer <A. 
Similarly 

C 0 = (2$ 0 ) im (2Vo)i n cos[mP*T 0 J= (24> 0 ) im (2F 0 ) in 

lim] [in] 

Xcos[w0*P— mP\T— Po)]= Y Y (— 1)« 1+32 

91=0 9 2=0 

XC m i qi C\ qi {P^) Ml R m - i HQ^) Ut Z n ~ Mi 

(i(m-l)] [j(n-l)] 

+ E E (-1) ^c m 2qi+1 c^(P^+ l 

91=0 9 2=0 

X ^- 291 - 1 (^ 2 ) 29 ^ 1 ^^ 292 - 1 . ( 12 ) 

Hence So and Co are two polynomial integrals of the 
system (10), but they are not independent. In fact 

5o 2 +Co 2 =(24> 0 )-(2Fo) n . (13) 

These integrals have terms of degrees a , b , c , d , in £, 
P, 2, Z, respectively, where a+&=w, c+d=w, and 
a-\-c= odd in So, while a+c=even in Co. It is seen that 
n is always even, because 7+6 and Ty are always even 
(paper A, pp. 279-280). 

It can be easily proven now that every polynomial 
integral of the system (10) is a polynomial of 4> 0 , Vo, 
So and C 0 . In fact, such an integral can be written as 
a polynomial of the quantities (6). If we express the 
powers of sines and cosines by means of sines and 
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cosines of multiples of x=P*(T— T 0 ) and y=Q*T we 
get terms of the form 

cos 

q( 2$ 0 )* {a+fi) (2Vo)^ +6) . ( u-2p 1 )xcos(l3-2p 2 )x 

sin 

cos 

X (y~2pz)y cos(8—2p A )y, (14) 
sin 

where we have cosine or sine whenever a or y is even 
or odd respectively; a, /?, 7, 8 are positive integers and 
pi, Ph P* are integers absolutely equal or smaller 
than a , 13, y, 8. The terms (14) are sums of terms of the 
form 

cos 

q'( 24>o)* (a+ * ) (2F 0 )* ( ^ a) (AiX±T jy ), (15) 

sin 

where 

A ,= (a— 2pi)±(fi— 2p 2 ), 

Tj=(y— 2p 3 )±(8— 2p t ), 

hence 

A &± Yjy = (A iP*± T&i) T- A iP*T 0 . 

As the integral is independent of T, all terms of the 
form (15) will be zero except those that have 

AiP*±T@*=0 9 i.e., Ai=±km, Fj—dakn. Then a+/3 
=km+ 2vi*, y+d=kn+ 2 v 2 *(k> 0 , vi*>0, v 2 *>0). 

Therefore the integral is a sum of terms of the form 

±g'( 2$o) i(2 ' 1 * ) (27 0 )* (2 '** ) 

[ COS “I 

(2*o) |fa "(2Fo)»*" (kmP h T 0 ) (16) 
sin J 

or ±^'(2^ > o)' ,1 *(2To) , ' 2 *X (polynomial in So, C 0 of degree 
k), i.e., it is a polynomial in 3> 0 , Vo, So, Co. 

If now in formula (7), ($„fi) includes a term £ a RPz y Z s 
with a-\-/3=m+ 2v h y+8=^n+2v 2 (u 1 >0, v 2 >0), then 
$y+ 1 contains secular terms of the form 

cos 

dbg(2$o) i(a+ ^(2T 0 )^^ 5) (■ mP'T 0 )T ; 

sin 

in fact if a = even, y=even we have terms of the form 
(paper A, p. 279) 

2^(24>o) i(ac+/5) (2Fo) i(7+5) J cosAiX cosTjydT 

= j(2* o)*<*W ) (27 0 )*<rH> 

X j " [cos(^4 ».r+ T ; *y)+cos(^4 {X— T jy)~]dT. 

Ai, Tj take all the values from (a+/3) to — (a+/3) and 
from (t+ 5) to —(7+5), respectively. Therefore, there 


are values of A i} Fj for which Ai= dt m, T,= ±n, i.e., we 
have secular terms of the form 

?(24>o)* ( ^2To) i(T+5) cos ( wPTo)r. 

If a=odd, 7= odd, we have terms of the form 

2</(24> 0 )R" + ^ ) (2F 0 )* (r+a) j' sin A x x sin FjydT 

= q( 24>o) i(a+ ^(2F 0 ) i( ^ 5) 

X y [cos(^4 iX— T 3 y) — cos(A X x+ Fjy^dT, 

hence we have secular terms of the form 

±q( 2<t> 0 )*<^ ) (2Fo)* ( ' y+5) cos (mP*T 0 )T. (17) 

If a = even, 7= odd we have 

2q(2$o)* (a+fi) (2Vo)* (y ~ hS) J cos A x x sin FjydT 

= q( 2<t> 0 )* ( ^ ) (2F 0 )* ( ^ a) 

X J [sin(r i y+i4^)+sin(r i y— Aix)~\dT, 

i.e., the secular terms are of the form 

±g(24> 0 )i ( ^K2Fo) i( ^ a) sin(wPir 0 )r, (18) 

and if a = odd, 7 = even we have 

29(2<fo) 1<<rf ' ,) (2Fo) J< - >+i) J sinAiX cosT, ydt 

= 9 (2<t> 0 )*<“+«(2Fo) 1(r+i) 

X J [sin(^ <a:+ r,-)»)+sm (^ &— Tjy)~\dT, 

i.e., we have secular terms of the form (18). 

Therefore if <*+7= even we have secular terms of the 
form (17) and if a+y = odd we have secular terms of the 
form (18) 

If a-\-f3=2m-\-2vi , y-\-d=2n-\-2v 2 (vi>0,v 2 >0) we 
also have secular terms of the form 

cos 

dz^(24>) i(a+ ^(2Fo) i(rfa) (2 mP'To)T, (18a) 

sin 

etc. 

If we calculate successively the terms 4> 0 , <t>i, 4> 2 * • - of 
the third integral, where <$>„ is of degree v-\~2 in £, R, 
z, Z, then the first secular term is in <t>„ with v=m-\~n— 2 • 
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If we replace <£ 0 by So or Co, where a+b = m } we get 
the integral 

S=S 0 +bS 1 +b 2 S 2 +- •, or C=C 0 +bC 1 +b 2 C i +- - . 
Then there will appear two secular terms in S 2 or C 2 : 
cos 

qi(2$ 0 ) * (-+ 2 ) ( 2 Vo) (mP'To) T 

sin 

cos 

+^(2^ 0 ) im (2F 0 )^ w+2 > b (mP'To)T, (19) 
sin 

unless it happens that some coefficients q are zero. 

In the terms (19) we have cosines for S 2 and sines for 
C 2 because in S 0 it is a+c=odd, hence in S 2 , S 4 - • • : 
a+c=odd, but in Si, S 3 --*: a+c=even, therefore in 
(Si,Fi): a+Y=even (Paper A, p. 278). Similarly in 
C 2 , C 4 * • • we have a+c=even, and in Ci, C 3 * • • : 
a+c=odd, hence in (Ci,Fi): o :+7 = odd. 

In order to eliminate the secular terms (19) we may 
take an integral of the form 

<P= <po~\-b<pi-\-b 2 <p 2 -\- • • *, (20) 

where 

cpo=x 1 (2^o) 2 +x 2 (2^o)(2Vo)+x z (2Vo) 2 , (21) 

and in the term <p m +n- 4 either S 0 , or Co is included; in 
fact <p k is of degree £+4 in £, R, z, Z. Further, as 
<Po 7 <p 2 , **• have a+c=even and <p u <p 3 , ••• have 
a+c=odd, S 0 is used in <p m+n _ 4 if w+w=odd, and Co 
if w+w= even. 

Then we define Xi, x 2 , x 3 so that the secular terms in 
<Pm+n- 2 disappear. In fact the terms (21) give also two 
secular terms of the form (19) in <pm+ n - 2 , which should 
cancel the corresponding secular terms due to So or Co. 
Thus we have two linear equations for the three un- 
knowns X\, x 2) x 3 . 

If £ 10 , # 20 , #30 is one solution of these equations then 
the other solutions are #i=#io+#i', # 2 = # 20 + # 2 ', 
#3=#3o+#3 / , where Xi , x 2 \ x 3 satisfy two linear homo- 
geneous equations. Then x\=*\x 2 =xi\ in fact we see 
that if the quantity 

^ , (2Fo) 2 =x 1 '[(2^>o) 2 +2(24>o)(2Fo)+(27o) 2 ] 

is included in <p 0 , this gives terms up to <i> 2 , but not 
higher order terms ; its contribution to the secular terms 
is zero. This means that if we set # 2 '= 2 * 1 ', x^=xi, we 
have a solution of the above system of two homogeneous 
equations in #/, x 2 , x 3 . But this system has only one 
solution with respect to x 2 /x\, x 3 fxi\ i.e., this is the 
solution given above. 

The secular terms of S 4 or C 4 are eliminated in a 
similar way by adding 

Qyi(24>o) 3 +y2(24> 0 ) 2 (2Fo)+y 3 (2$o)(2 F 0 ) 2 +y 4 (2F 0 ) 3 ] 


in (p 2 . Then both this quantity and So or Co give three 
secular terms in (p m+n of the form 

cos 

±?(2$o)* (m+4) (2F 0 )* ,n < (wF^o)F, 
sin 

cos 

d= g(23> 0 ) * (,w+2) (2 F 0 ) * (n+2 > (mP'To)T, (22) 

sin 

cos 

±®( 24> 0 )*’ n (2F 0 )* ( "+ 4) b (tnP*To)T. 

sin 

If we require that these two sets of secular terms 
cancel each other, we have three equations for the four 
quantities y 1} y 2 , y 3 , y 4 . It can be proved as above that 
if y 10 , 720 , ^ 30 , 740 is one solution of these equations, 
then the general solution is yi=yio+yi, y2=y2o+y2, 
y*=yzo+yz, y*=y4o+y4, where yi'=\y2=\yz=yi- 
In general v? m+n+ 2 fc-4 includes S 2 k or C 2 k, that has 
k+1 secular terms 

cos 

±q(23> 0 ) Hm+2k) (2Vo) in . (mP*To)T, 
sin 

(23) 

cos 

=t (24> 0 ) * m (2 V 0 ) H w+ 2 *o . (mP*To)T 
sin 

and these are eliminated by adding the quantity 

«i(2$o)* + H“«i(2$o)*(2Fo)+ • • • -\-Wk+ 2 ( 2 Vo) k+l 

in cp 2 k- 2 , which gives again £+1 secular terms in 
(Pm+n+2k-A) by setting the coefficients of the secular 
terms of <pm+n+ 2 fc -4 equal to zero, we have k+\ equations 
for the k-\-2 unknowns coi, co 2 , • • -co* + 2 . If o>io, co 2 o, 
• * -o>k+ 2,0 is one solution, then all other solutions are 

0>l = COiO+w/, CJ 2 = C0 2 o+W 2 / ,* * •COJfc-|-2 = WA:+2,0^ _ ^A:4-2 , , 

where 

co 2 co 3 / 

Oi\ = = = * ’ * ~&k+2 • 

C^\ C k+ 1 2 

In order to eliminate the secular terms of the form 
(18a) in $ 2m +2n+2k-4 we add a quantity 

So 

[aj 1 (24 > 0 ) A:+1 +aj2(24> 0 ) A: (2F 0 )H — • +coa : + 2 (2F 0 )* +1 ] 

Co 

in 4> m+n+2 fc-2; again we take S 0 if m+n= odd, and Co is 
w+^=even. Then both this quantity and So (or Co) 
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give k + 1 secular terms of the form (18a) in $ 2 m+ 2 n+ 2 k- 4 : 

cos 

d= <?(2<l>o) * i2m+2k) (2 F 0 ) * (2n) . (2 mP*To)T, 
sin 


In the same way the secular terms of the most general 
form 

cos 

±q(2$ 0 ) i{a+fi) (2Vo)^ ) . ( jmP>T 0 )T 
sin 


±q(2<i>o) H2m+2k - 2) (2Vo) H2n+2) 


cos 


(2 mP*To)T, 
sin 


(24) 


cos 

±g(24> 0 )i (2w) (2Fo) i(2n+2A) (2 tnP*T 0 )Ti 

sin 

by requiring that these terms cancel each other we get 
k+1 equations for the k-\-2 unknowns «!,•••, <o*+ 2 . 
The general solution is again of the form 

Wl=a>lO+Wi', * * * 03b+2= W*+2,0+^+2 / 

where duo, • • • , uic+ 2,0 is a special solution, and 
=u 2 ' /C k+l i~ • • • =<0*4.2'. 

In fact the integral 


So 

7o'= (2F 0 ) k+1 

Co 


does not give any higher-order secular terms beyond 
those given by So or Co; in particular it does not give 
any secular terms of. the form (24) in $ 2 m+ 2 n+ 2 k- 4 - If I\ 
is the next higher order term after 7 0 ', then 


(I 1 ',F 0 )=-(Io',F l )=- 



*+1 


S 0 

-2(A+l)(Po,Pi) (2F 0 )*, 

Co 


s s 

7/= \2Fo) k ^+2(k+\)F l °(2P 0 )*, 
Ci Co 

because 

- 1 (F 0 ,F 1 )dT=F 1 . 

Similarly 


are eliminated by adding terms of the form 


So^ 1 

[O 1 (2* 0 )* l +O 2 (2*o)*(27o)+ • • • +U k+2 (2Vo)^ 

Co*” 1 

m (p (J—l) (m+n)+2k—2‘ 

Therefore we can always eliminate the secular terms 
in the third integral, written in the form (20). 

Now we may distinguish a number of cases. 

(a) If m+n> 4, then <3> 0 depends only on <t> 0 and V 0 ; 
So (or Co) has as a coefficient some positive power of b , 
therefore if b is small the influence of this term is small. 
This means that if m-\-n > 4 the form of the third 
integral depends mainly on <f> 0 and Vq=Fq—$o. 

(0) If m+n= 4, then 4> 0 includes both £i(24> 0 ) 2 
+o: 2 (2 <I>o)( 2 Fo)+^ 3(2 Fo ) 2 and Co. In fact, then the 
secular term of <f> appears first in <t> 2 . Then m = n= 2, 
i.e, 7>*=C*. 

(7) If m+n< 4, we write <p in the form 

^ = S 0 + b\jXi(2$o) 2 + x 2 (24> 0 ) (2 F 0 ) + x 3 (2 F 0 ) 2 J 

+^(... )+.... (25) 

Then the secular term of 4> appears in 4>i, i.e., 
4Q-P=0, or P*/Q>=2/1 , m= 1, n=2. 

It is seen that m-\-n cannot become smaller than 3, 
as 4>i already is of the third degree. In fact, n is always 
even, but neither m nor n becomes zero. 

This distinction between the cases a, 0, and 7 shows 
that higher order commensurabilities (i.e., if * s 

not equal to 1/1 or 2/1) do not give any appreciable 
change of the over-all form of the orbits from the general 
irrational case, if b is small, while if P*/Q*= 1/1 or 2/1 
we should expect to find important differences in the 
general forms of the orbits. This is, in fact, what has 
been observed in the calculated orbits. 

2. GENERAL RATIONAL CASE 

Some orbits have been calculated in the case (I) 
pi/Qi= 4/1 (m+n= 5) and (II) P*/Q*= 2/3 (m+n= 5). 
In these cases $0 gives secular terms in <£ 3 . 


6" S 

U = 2 (2F 0 ) k+1 +2(k+ l)Fi 1 (2Fo) k } 

C 2 C x 

because F 2 =0, and finally 

S m+n—2 Rm+n — 3 

I m +n-2= , (2Fo) k+l +2(k+l)F l (2Fo) k 

^ m-fn— 2 C m 4-n— 3 

and this term is of degree 2w+2w+2£ but has no secular 
terms of the form (24). 


Case (I) P= 16Q 

In this case w= 1, n= 4, m+n= 5, therefore we must 
use the integral So, or 

S 0 = —So/4Q*=Q 2 t;z 4 — 6Q£z 2 Z 2 

+ & 4 -RzZJ+QRz*Z (26) 

of the unperturbed system. 

In the general form of the expansion of the third in- 
tegral in powers of b , some third-order terms (Barbanis 
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1962) have zero divisors. If we multiply $3 (the coeffi- 
cient of f> 3 )_by (P—160/180 and then set P=160 we 
find again So- This should be expected because the only 
integral of fifth degree of the system (10) in this case 
isVS’o or a multiple of it. 

If then we calculate Si by the formula 

(<Si,P 0 )+(<§o,Pi) = 0, (27) 

we get 


. (P 2 +P£ 2 ) 

Si = (- 1 1Z 4 - 540 2 Z 2 + 290z 4 ) 

2560 


(P 2 -P0 

+ 

7680 

2£RzZ 

+ (Z 2 — 20 2 )— 0 2 z 4 — — 

3 Q 6 


[8Z 4 - 1440 2 Z 2 + 160z 4 ] 


0 « 


(Z 2 -£>z 2 ) 

[Z 4 +40 2 Z 2 +0z 4 ] , 

240 


(28) 


and in the calculation of S 2 there appear the secular 
terms 


A UP 
1920 1920’ 

where 

A = T(2<t>o) 3,2 (2F 0 ) 2 cosP 1,2 P 0 

and 


(29) 


P= P(2<I>o) 1/2 (2Fo) 3 cos P^To. 


On the other hand, the terms 

x 1 (2^o) 2 +x 1 (2^o)(2Vo)+x 3 (2V 0 ) 2 (30) 

in ipo give 

8^id>o4>i+ 4x 2 ('t ) o F id- V od>i)+ 8.T3 V oV 1 (31) 

in <p h where Fi=Fi— $1, and 

4^Zi(4 > x 2 — f- 2 $o4 > 2)+4x2(4>o F 2 _ t _ 4’i F id - I^F 0) 

+4x 3 (Fi 2 +2F 0 F2) (32) 

in ip2, where F 2 = —$2 (because V=F— <f> and P2=0). 
In <p 3 the terms (30) give the secular terms 

xi A x 2 / A B \ x 3 B 

+-( +— ) . (33) 

96<2 4 (A 192 192/ 96<2 4 

Adding the secular terms (29) and (33) we must have 
zero, i.e., 

A B 

(Q+ 2 x 1 -x 2 )+ (- 11Q+x 2 - 2xt) = 0. (34) 

192Q 4 192 Q 4 

Hence one solution for x h X 2 , x% is 

*io =— 10, *2o=0, *3o= — (11/2)0 (35) 

and the other solutions are *i=*io+*i / , *2=*2o+*2 / , 
*3= *30+ *3', where *i'=£*2 , = *3 / , in agreement with 
the general theory. 
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Table I. Characteristics of the orbits represented in Fig. 1. 



£o 

zo 

*0 

Zo 

2 F 

in 

(2=0) 

£max 
(2 = 0) 

& 

Zi 

£2 

z 2 

a 

0. 

0. 

-0.09850 

0.07482 

0.015300282 

-0.079 

0.074 

-0.070 

±0.242 

0.085 

±0.234 

b 

0. 

0. 

0.09850 

0.07482 

0.015300282 

-0.079 

0.075 

-0.069 

±0.242 

0.085 

±0.234 

c 

0. 

0. 

0.05120 

0.11260 

0.015300199 

-0.044 

0.036 

-0.022 

±0.365 

0.057 

±0.358 

d 

0. 

0. 

-0.05120 

0.11260 

0.015300199 

-0.044 

0.036 

-0.022 

±0.365 

0.058 

±0.358 


Therefore the third integral is 

*>= -$G(2* 0 ) 2 - (11/2)0(2 F 0 ) 2 

+b(-4Q^i-4AQV 0 Vi+S 0 ) 

+bK-2Q(^+2^2) 

-22Q(V 1 *+2VoV 2 )+S 1 )+ • • - (36) 

and no secular term appears now in <p 3 . 

Four orbits were calculated in the potential field (1) 
with P= 1.6 (10 7 yr)" 2 , Q= 0.1 (10 7 yr) -2 , b= 0.2 kpc -1 
(10 7 yr) -2 ; C is not needed in our calculation. 

This field does not represent any real stellar system; 
in fact it should represent approximately a prolate 
stellar system. The results, however, are typical of any 
case of higher-order commensurability between P* 
and QK 

The four orbits are represented in Fig. 1 (a, b, c, d). 
Their characteristics are given in Table I: In this table 
£o, z 0 , Ro, Z 0 represent the initial coordinates and 
velocities; £ ra in(z=0), £ ma x(z=0) are the minimum and 
maximum values of £ when the curve crosses the £ axis; 
(£i,zi) and (£ 2 , 22 ) are the angular points of the boundary; 
and F is the value of the energy. The units used are 
kpc for £ and z, kpc (10 7 yr) -1 for R and Z, and kpc 2 
(10 7 yr) -2 for F. We have 

2F— R 2 + 16Q£ 2 +Z 2 +Qz 2 - 26£z 2 = z> 0 2 

and the measure of the initial velocity is approximately 
the same: v 0 = 0.1 23694 kpc (10 7 yr)'" 1 . 

The calculations of the orbits were made by means of 
Runge-Kutta method, in double precision. Whenever 
a step 0.02 (10 7 yr) was used in integration, the accuracy 
of the energy integral was within 1 unit in the ninth 
decimal place for 3X10 9 yr (15000 steps). Whenever 
the step was 0.05 (10 7 yr) the energy integral was con- 
served within two units in the eighth decimal place for 
3X10 9 yr. A comparison of the individual values of 
coordinates and velocities for the same set of data and 
different steps showed that at least five decimals were 
the same for an interval of 2X10 9 yr. 

When single precision with step 0.05 (10 7 yr) was used, 
the error in the energy constant was in one case one unit 
at the sixth decimal place after 1.1X10 9 yr. In another 
case, when a step 0.1 (10 7 yr) was used, one unit at the 
sixth decimal place was lost after 1.7X10 9 yr. 

A comparison of the individual values of coordinates 
and velocities between single and double precision 
showed an agreement within two units in the fifth 
decimal place for 2X10 9 yr, when a step 0.1 (10 7 vr) 


was used, and within one unit in the fifth decimal place 
when a step 0.05 (10 7 yr) was used, although in this 
case the energy integral was conserved with less 
accuracy. 

From Fig. 1 and Table I it is seen that the boundaries 
of the orbits are curvilinear “parallelograms” whose 
right and left boundaries are slightly concave to the 
right. The apices of these curvilinear parallelograms 
were found from an inspection of the tables of the values 
of the coordinates for every 0.05 X10 7 yr, which are not 
given here. These apices lie, within an accuracy of 
±0.001, on the curve of zero velocity, i.e., the curve 
found by setting R==Z= 0 in Eq. (2). In our case the 
curve of zero velocity 

Pi; 2 +Qz 2 -2b&=vo 2 (37) 

is represented by a dashed line in Fig. 1. Its deviation 
from the ellipse 

P?+Qz 2 =Vo 2 

is measured by the maximum value of 26£z 2 /z> 0 2 . If we 
write 

co =2 ££z 2 , 

then 

2&£(i> 0 2 -P£ 2 ) 

Q- m ’ 

and this is maximum if 

46P£ 3 -3P(7£ 2 +ti 0 2 (7=0. 

For P= 0.4, (7=(U, 6=0.2, and i> 0 2 =0.0153, the 
maximum value of gj is co mnx = 0.0030. Therefore the 
maximum value of co/» 0 2 is 20%. 

The curve of zero velocity intersects the £ axis at the 
points £= ±z>o /P i= ±0.098, and the z axis at the points 
z=±VC*= ±0.391. The maximum value of z is for 
£ given by the equation 

bP?-PQZ+bv 0 2 = 0. 

The solution of this equation that is smaller than 
vo/ p* is 

t 2bv 0 2 /PQ bv 0 2 ^ ( b 2 v 0 2 ^ 

{1+[1-(46W/P0 2 )] ! } PQ\PQ> + / 

For P=1.6, (7=0.1, 6=0.2, do 2= 0.0153 we have 
2max= 0.399 for £=0.020. 
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Table II. Approximate values of the third integral. 



Initial 

<pQ 

Max 

Min 

Initial 

<po~i~b<pi 

Max 

Min 

Initial 

T" b -f~ bP 
Max Min 

T 

a 

-0.22 

-0.19 

-0.24 

-0.219 

-0.217 

-0.219 

-0.2188 

-0.2185 

-0.2191 

4X 10 s yr 

b 

-0.22 

-0.20 

-0.25 

-0.219 

-0.217 

-0.219 

-0.2188 

-0.2186 

-0.2191 

5X10 9 yr 

c 

-0.89 

-0.87 

-0.99 

-0.888 

-0.884 

-0.956 

-0.8894 

-0.8858 

-0.8938 

3X 10 9 yr 

d 

-0.89 

-0.80 

-0.99 

-0.888 

-0.884 

-0.957 

-0.8894 

-0.8846 

-0.8931 

SX10 9 yr 


Equation (37) gives 

z*=(v 0 *-m/(Q-2bO. 

In order that the curve of zero velocity should be 
closed, so that the orbits should not extend to infinity, 
we must have Q— 2b £>0 whenever v 0 2 — P£ 2 >0. Hence 

Q-2bv 0 /P*>0, 

or 

b<QP'/2v 0 . 

In the above case b <0.51. 

If 5=0, the boundaries are exactly parallelograms 
with sides 

£= ±Po/ 4 Q\ z= ±Z 0 / Q *, 
i.e. 


orbits a, b: £=±0.078, z= ±0.237; 
orbits c, d: £= ±0.040, z= ±0.356. 

Therefore the deviations of the actual boundaries 
from the corresponding parallelograms of the unper- 
turbed cases are rather small. 

The third integral can be written 

<po+5<pi+5V 2 H — •, (38) 

where 

<Po= -^(i? 2 +16^ 2 ) 2 -(ll/2)G(Z 2 +Q 2 2 ) 2 , (39) 

<Pi = Q 2 & 4 - 6(9£s 2 Z 2 + £Z 4 — RzZ z +QRz z Z 

+l(R 2 +16Q?)(7Qiiz 2 -ZZ 2 +RzZ) 

- (ll/3)(Z 2 +Qz 2 )(Q& 2 - $Z 2 +RzZ) (40) 

and 


(. R 2 + 16(?£ 2 ) (R 2 - 16<2£ 2 ) 2 £RzZ 

1 — (— 11Z 4 — 54Qz 2 Z 2 +29Q 2 z 4 )+ (Z 4 -18Q2 2 Z 2 +2<2 2 2 4 )+ (Z 2 -2Q 2 2 ) 

256Q 2 96Q 2 3Q 

_0« 2 4 1 — —(Z 4 +4Qz 2 Z 2 +Q 2 z 4 ) [(7Q& 2 — £Z 2 +RzZ) 2 + U(Q(z 2 -iZ 2 +RzZ) 2 ] 

6 24Q 2 18() 


- J[7? 2 + 16 Q ?- 1 1 (Z 2 +Qz 2 )] ■ 


'z* 11 R 2 Z 2 9 | 2 Z 2 2lR 2 z 2 | 2 z 2 2 £RzZ\ 


+ 


+ 


, 2 160<2 2 10 Q 160(2 10 5(2 


-)■ 


(41) 


The values of <po, 5<pi, v?o+5<pi, 5 2 <p 2 , <po+5<pi+5 2 v? 2 
were printed every 10 8 yr. Table II gives the initial 
values of <po> ^o+5<^i, ipQ-\~b(p\-\-b 2 (p 2 and their maxi- 
mum and minimum values for the time intervals given 
in the last column, in units of 10 -4 kpc 2 (10 7 yr) -4 . It is 
seen that if we include terms up to <p 2 in the third 
integral its accuracy is fairlv good of intervals of at 
least some billion years. 

Further a great similarity between orbits a and b and 
between orbits c and d is to be remarked (Tables I 
and II and Fig. 1). This indicates that two orbits 
with velocities symmetric with respect to the axes 
£ and 2 have the same general behavior (the same 
boundaries, etc.). 

Case (II) 9P=4£) 

Five orbits were calculated in this case with P=0A 


(10 7 yr) -2 , (9=0.9 (10 7 yr) -2 and 

orbit a: 5=0.2 kpc -1 (10 7 yr) -2 , £o=zo=0, 

Po=0. 11260 kpc (10 7 yr) -1 , Z 0 =0.05120 kpc (10 7 yr) -1 ; 
orbit b: 5=0.4 kpc -1 (10 7 yr) -2 , £o=zo=0, 

Po=0.11260 kpc (10 7 yr) -1 , Z 0 =0.05120 kpc (10 7 yr) -1 ; 
orbit c: 5 = 0.4 kpc -1 (10 7 yr) -2 , £o=Zo=0, 

Po= 0.05 120 kpc (10 7 yr) -1 , Z 0 =0.11260 kpc (10 7 yr) -1 ; 
orbit d: 5 = 1 kpc -1 (10 7 yr) -2 , £o=zo=0, 

Po=0.11260 kpc (10 7 yr) -1 , Z 0 =0.05120 kpc (10 7 yr) -1 ; 
orbit e: 5=2 kpc -1 (10 7 yr) -2 , £o=zo=0, 

Po= 0.1 1260 kpc (10 7 yr) -1 , Z 0 = 0.05 120 kpc (10 7 yr)" 1 . 

In all these cases the velocity is the same : Vo = 0. 1 23694 
kpc (10 7 yr) -1 . 
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Fig. 2. Orbits in the case 9P = AQ. 


The calculations were made again by means of 
Runge-Kutta method in double precision. The energy 
integral 2F( = 0.015300199) was conserved with the 
following accuracy for the corresponding time intervals : 


orbit a: 

3X10" 9 

for 

3X10 9 yr 

orbit b: 

16X10- 9 

for 

12X10 9 yr 

orbit c: 36X10 -9 

for 

9X10 9 yr 

orbit d: 

36X10- 9 

for 

9X10 9 yr 

orbit e: 

9X10- 9 

for 

5X10 9 yr 


The forms of the orbits are represented in Fig. 2(a-e). 

Orbit a (6=0.2) remains quite near the periodic orbit 
of the unperturbed case (6=0) for a few billion years; 
but as the perturbation increases the orbits deviate 
more and more from the unperturbed case; orbit e does 
not resemble the periodic case at all. 

In all cases, the orbits do not fill all the space inside 
the curve of zero velocity 

0.4£ 2 +0.9z 2 —26£z 2 = 0.0153, 

which is represented by a dashed line in Fig. 2. In order 
that this curve should be closed we must have 6 <2.30; 
case e is near the limiting case. 


In cases a, b, c, the boundary of the orbit is approxi- 
mately a parallelogram with sides 

Z=±Rq/P\ z=±Z 0 /Q\ 
i.e., 

orbits a, b: £=±0.178, z=±0.054; 
orbit c: £=±0.081, z= ±0.119. 

The actual boundaries are curved in the same way as 
in cases c and d but in a smaller degree. The character- 
istics of the orbits are given in Table III. The boundary 
is symmetric with respect to the £ axis. Both left and 
right side of the boundary are concave to the right. The 
points where these sides intersect the £ axis are given by 
£min(2=0) and £ max (z=0). The upper and lower sides 
join the left and right sides at the usual angular points 
of the boundary (£i,zi) and (£ 2 ,z 2 ). These points lie on 
the curve of zero velocity with an accuracy of about 
±0.001 in both coordinates. 

The points (£ 3 , 23 ) and (£ 4 , 24 ) are the minimum and 
maximum of the upper boundary, and correspondingly 
the maximum and minimum of the lower boundary. The 
accuracy of the values £3 and £4 is only about ±0.01. 

The point (£4,24) is always a regular point of the 
boundary. The same seems to be the case for the point 


Table III. Characteristics of the orbits represented in Fig. 2. 


£min £max 



( 0 = 0 ) 

( 0 = 0 ) 


2 l 

£2 

Z2 

£3 

Z3 

£4 

24 

b 

- 0.179 

+ 0.177 

- 0.172 

± 0.056 

± 0.181 

± 0.054 

- 0.10 

± 0.050 

± 0.12 

± 0.060 

c 

- 0.079 

+ 0.089 

- 0.073 

± 0.116 

± 0.092 

± 0.119 





d 

- 0.178 

+ 0.172 

- 0.158 

± 0.067 

± 0.184 

± 0.058 

- 0.11 

± 0.045 

± 0.08 

± 0.074 

e 

- 0.174 

+ 0.166 

- 0.104 

± 0.091 

± 0.190 

± 0.081 

- 0.08 

± 0.055 

± 0.09 

± 0.110 


orbit a. £min — 0.178, £ max — I 0 . 179; 2 max and 2min vary from zfc0.052 to :i0.055. 
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($ 3 , 23 ) in cases b, c, d; in case e, however, this minimum 
seems to be an angular point of the boundary. 

The curves of zero velocity intersect the £ axis at the 
points £= ± v 0 /P*= ±0.196, and the 2 axis at the points 
z=±i>o/(? i = ±0.130. The maximum and minimum 
values of z are: 

case a (6=0.2): z=±0.131 for £=0.009; 

cases b, c (6=0.4): z=±0.131 for £=0.017; 
cased (6=1): z=±0.134 for £=0.045; 

case e (6=2): z=±0.151 for £=0.114. 

The deviation of the curve of zero velocity from the 
ellipse 

0.4£ 2 +0.9z 2 =tt 0 2 

is again given by the maximum value of 2 6£z 2 /fl 0 2 , which 
is: case a: 3.5%, cases b, c: 7.5%, case d: 23%, case e: 
82%. It is seen that in the last case the deviation is 
considerable; however, even in this case the existence 
of an accurate boundary of the orbit indicates that a 
third integral exists, if not for all times, at least for 
times of the order of many billion years. 

The exact form of the third integral was not calcu- 
lated in this case. 

3. the case P/<3=4/l 

In this case the third integral <p can be given by 
Eq. (25), where S 0 is replaced by 

So=So/2Q*= £(Qz 2 — Z 2 ) + RzZ. (42) 

This is the limit_of $i(4(7— P)/2 when P — > 4 Q. 

If we calculate Si by the formula (27) we find 

1/z 4 £ 2 Z 2 z 2 R 2 2 £PzZ\ 

&=— ( £ 2 z 2 + ). (43) 

2\2 Q Q Q / 

This is the limit of 

(4 Q-P) (2 P+Q) 


4> 2 


iPQ(P-Q) 


(2<to)(2Fo) 


when P — » 4(). 

Now $ 2 has the secular terms 


where 

and 


-5^/64<2 3 +5P/16() 3 , 

A = T( 2$o) 3/2 (2To) cosP 1/2 Po 
B = T (24 > 0 ) 1 12 (2 V o) 2 cosP^Po. 


(44) 


On the other hand, the term £i(24> 0 ) 2 +£ 2 (2$o)(2Fo) 
+xz(2V 0 ) 2 of <pi gives the secular terms 


X\ A /—A B\ x%B 

h^2i 1 ) 

Q V 2 Q 2 Q/ Q 


( 45 ) 


Adding the terms (44) and (45) we must have zero. 
Hence 


-( — 5 + (AQ 2 x\— 32 Q 2 x 2 ) 


64 Q 


B 


16 Q 3 


-(5+80 2 .v 2 -16(7 2 x 3 ) = O. (46) 


One solution of this equation is 

1 3 7 

*10= , *20= , *30= , 

64@ 2 16() 2 32Q 2 

and all other solutions are *i=*io+*/, *2= *20+ *2', 
*3=*3o+*3', where */= \xj=xi. 

Therefore the third integral becomes 

6/z 4 £ 2 Z 2 z 2 R 2 2£PzZ 

V — £((?z 2 — Z 2 ) + RzZ — ( £ 2 z 2 H 

2 \Q Q Q Q 


1 


32Q 2 


-[(P 2 +P£ 2 ) 2 + 12(P 2 ±P£ 2 )(Z 2 +()z 2 ) 


— 14(Z 2 ±@z 2 ) 2 ] H . (47) 


This can be written 

( £ 4 £ 2 z 2 z 4 £ 2 P 2 £ 2 Z 2 

(P=Q£z 2 -£Z 2 +PzZ+6 

l 4 4 32 SQ 40 


5i? 2 z 2 

-1 1 

7z 2 Z 2 

%RzZ 

R* 

3i? 2 Z 2 

, 7z ‘l 

1 l 

16(3 

160 

Q 

640 2 

160 2 

320 2 l 


64 Q 2 


+ • 


-(Po 4 +12Po 2 Z 0 2 -14Zo 4 )+- ■ •, (48) 


where P 0 , Z 0 are the components of the initial velocity 
i’o, when £o=zo=0. 

The same result can be found also in the following 
way. The coefficient 4> 3 of 6 3 (Barbanis 1962) can be 
written after some operations: 

Q r 96(23 > o) 504(2 F 0 )l 

4>3 = So 

(4 Q-P) 2 (P-Q)P L(9P-40) (P— 160 J 

1 

+ terms with coefficient- 


(4 Q-P) 

If we multiply the integral by J(4 Q—P) the only 
critical terms in |[03(4(9— P)] are 

QS 0 r 96 ( 2 ^o) 504(2 F 0 )“| 

. (49) 

±9P-40 P-160J 


2(4Q-P)(P-Q)PL9P-4Q P- 1 6Q . 


in cp 2 . 
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But if we add 


0 pyi(z<Por 

(p-0)pL9P-40 p- i6yJ 

in £[$2(40— P)], we get the further terms 

4 Q ryi$i(2<f> 0 ) y2V 1 (2V 0 )’ 


( P-Q)P\ 


pi^i^o) y^y i\2v ojn 

L 9P-40 p-i 6 yJ 


i n i[$3(40— P)], and these give the further critical 
terms 


40$i r yi(2w 0 ) y2(2Vo)~\ 

, (50) 

(P-0)Pl_9P-40 P-160J 

because Fi= — £z 2 — $i, and $i is critical. 


h If we add (49) and (50) and set a§ o /(40 — P) — > §0i 
when P — > 40, we must have zero. This is effected if 
y i= — 6 and y 2 = —63/2. 

Thus ¥?o=*§o and ip\ is the limit of 

(4 Q-P) Q r 6(2<t> 0 ) 2 63(2 F 0 ) 2 ~] 

$2 1 . 

2 (P— 0)P L9P— 40 2(P— 160 J 

By adding and subtracting [(2P+0)/4P0(P— 0] 
X(2$ 0 )(2F 0 ) we find for P=40 

3 (24> 0 ) 2 7(2 F 0 ) 2 

W = Si (24> 0 ) (2 Vo) + , 

160 640 320 

i.e., the same result as above. 

In the same way as above we can eliminate the secular 
terms from <p 3 , etc. The general form of <p 2 is 


r<f> 3 (4 0~P) 4 0 

<p 2 = lim 1 

p-4qL 2 (P— 0P 


64> 1 (24> 0 ) 63Fi(2F 0 ) 


9P—4Q 2(P— 160 


1 [13 75 55 

+/(4 ) o,Fo,>So) = 1 — 0£ 3 z 2 H — £ 3 Z 2 £ 2 RzZ 

24 0 1 4 4 4 


61 85 85 35 25 5 25 1 

— €*v+ £P 2 Z 2 RhZ 0£z 4 — 21 £z 2 Z 2 H £Z 4 — Pz 3 Z PzZ 3 +/(4> 0 ,Fo,So), (51) 

16 160 160 2 20 2 20 I 


where /o is an arbitrary function of $o, F 0 and So. 

The secular terms in <p 3 due to the explicitly written 
terms in (51) are 

475P/61440 9 ' 2 

where P= P(24> 0 )(2F 0 ) 2 s^P^Pq. If we add 

/(^o,Fo,So) = ^o[xi(2<I>o)+X2(2Fo)] (52) 

in <p 2 , we find also the secular terms (x\—x 2 )B/8Q 3 ' 2 . 
Hence if 


X\ x 2 — — 


475 
76 80 3 


the secular terms are eliminated. 

We take £io=0, x 2O =475/7680 3 and any other solu- 
tion is Xi=Xi 0 +Xi, x 2 =x 2 q+x 2j where Xi=x 2 . If 
we set 

/(* 0 , Fo,So) = (475/768<3 3 )(So(2 V 0 ) 


(53) 


in formula (51), we find 

1 (13 75 55 61 

<P 2 -= — 0 3 s 2 +— ?Z 2 ?RzZ £PV 

240 2 1 4 4 4 16 

85 85 85 

H £P 2 Z 2 RhZ 0£z 4 - 2 1 £z 2 Z 2 

160 160 32 


75 395 75 

£Z 4 H QRz 3 Z-\ RzZ z 

320 32 320 


(54) 


The third integral permits the calculation of the 
boundary of the orbit on the £,z plane. The method 
used is the same as in paper A. 

The energy integral (2) together with the third inte- 
gral (48) define a surface ((p,F) in the four-dimensional 
phase space; the projection of this surface on the plane 
£,z defines the space occupied by the orbit on this plane. 
The equation of the boundary of the orbit is then found 
by eliminating R and Z between Eqs. (2), (48) and 


dip dip 

J= — Z P=0. 

dR dZ 


(55) 


We will derive presently the equation of the boundary 
in first approximation with respect to 6; in all the calcu- 
lations below the terms of degree higher than the first 
in b are always omitted. 

Equation (2) gives 

R 2 = — Z 2 — 40£ 2 — 0z 2 +26£z 2 +fl o 2 , (56) 

and if we set this value in Eq. (48) and write 

(l/320 2 )(Po 4 + 12P 0 2 Zo 2 - 14Z 0 4 ), (57) 

we find 

Rzi 


J±\ 

Z(Z*-Qz*)+b 

| 25Z 4 

?z 2 

9z 2 Z 2 

q) 

l 640 2 

~~2Q~ 

32Q 

5» 0 2 Z 2 

_J_ 

3£V 23z 4 
1 1 

1 lfl 0 2 z 2 

v 0 * 

12 


' 320 2 

T 1 

2 64 

j 

32 Q 

64<2 2 

2 



(58) 
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On the other hand, 


[ / £ 2 2 2 5R 1 5 Z 2 

J=z{Z 2 -R 2 )+2^RZ+b\ RZ( 


5 Z 2 \ 
! 40 2 / 


and if we insert in this equation the values (63) and (65) 
we find the equation of the boundary of the orbit on 
the £, z plane: 


l \4Q 4 Q 16<2 2 4 Q 2 

= 0, (59) 


£ 2 z 2 


£z 

+-(R*-Z 2 ) 

<2 


-((2 2 z 4 -4CV$ 2 - 2<2»o 2 z 2 +»o 4 ) 


and if we set the value (56) in this equation we find 
in first approximation 


(* 2 +z 2 ) 

+45£ 2 X 6 (4(9£ 2 +3<2z 2 -V) 

-?(4Qe+3Qz 2 -v<?) (.4Qe+Qz 2 ~Vo 2 ) 2 


+26£ 


< 


2?z( 2£+5 


25Z 2 

? 

9z 2 

5V 1 

160 2 

Q 

1 

16(2 

160 2 ) 


4(£ 2 +z 2 ) 2 4 

— 25z 2 (» 0 2 — 2Qg‘—Qz 2 ) 2? 9z 2 5V 


= z(-2Z 2 -4<2£ 2 -(?z 2 -Ho 2 ) 

+5$ 


16Q 2 (£ 2 +z 2 ) (2 8(2 8g 


2£ 2 9z J 5V| 

1 

0 8(2 80 2 J 


(2 Z 2 V\ 


£ 4 z 2 (40£ 2 +30z 2 -V) 2 1 

40£V . (67) 

>( H4£ 2 +3z 2 ) 

\Q Q> 

. (60) 

0(f-+z 2 ) 2 1 


If we eliminate RZ between (58) and (60) we find 

/ 75 5£\ / 

Z 4 ( )+Z 2 ( 2(^+2 

V 320V V 

57z 2 5 V 


320 

5? 


-z 2 ) 


V 2 


5£( 57z 2 5V\ 

+- -2? + )+2(2^z 2 +ez 4 - 

gl 16 8 QiJ 

i 87 z* v 0 2 z 2 do 4 1 

+5£ -4£V 1 1 12 =0. (61) 

I 32 <2 320 2 J 

If we set 5=0, we find the zero-order solution for Z 2 : 

z 2 (V— 2Q?— Qz 2 ) 


2 (0 ) 2 = 

and if we set 
in Eq. (61), we find 
Xb=- 


2(e+z 2 ) 

Z 2 =Z(0) 2 +5X(, 


(62) 

(63) 


E 

87 z 4 V 2 » 0 4 

Afc2„2_l 

2(£ 2 +z 2 )l 

t-c Z “1 

32 <2 32 Q 2 


Q- 


2(0, 2 


57z 2 5V\ 75Z(o) 4 

-)+ 


v 57z 2 5t> 0 2 \ 

- 2£ 2 + - 

V 16 8(22 


(64) 


I* ^ 

5 R 2 

5Z 2 

1 4<2 4(2 

1 

16<2 2 

~4Q 2 


-z 2 (R 2 


’-M'-f)' 


In zero-order approximation this equation gives the 
following solutions: 


(68) 

(69) 

(70) 


(71) 


(a) Z(o) 2 — 0, 

(b) £(o> 2 =(V-<2z 2 ) 2 /W, 

(c) £(o, = 0. 

(a) If Z(o, 2 =0, we set 

z 2 =bY b 

in (67) and find 

(V-4<2£ 2 )(VF 6 +4$ 2 X 6 )=0, (72) 

where now 

1 r »o 4 I 

X„=- +G . 

2£L 32<2>' 2 -1 

Hence if V— 40£ 2 5^O, 


(73) 


2r v -| 

= — +o = 

fn 2 !— 320 2 J 


5Zo 2 (3Z 0 2 -27?o 2 ) 


and 


VL 32<2 2 J 16 QW 

55Z 0 2 (3Zo 2 -2i?o 2 ) 


0 \ 16 8(2/ 32<2 2 

Then the relation (56) gives 

R 2 = -Z m 2 -4Q?-Qz 2 +v 0 2 +b(2& 2 -X b ). (65) 

From Eq. (59) we find 


160V 


-f. 


■« (74) 
(75) 


( 66 ) 


This equation represents a parabola symmetric with 
respect to the axis £ and going through the origin. If 
b > 0 and 3Z 0 2 — 2Z? 0 2 >0 the parabola is to the right of 
the origin; if 3Z 0 2 — 2R 0 2 <0 the parabola is to the left of 
the origin; and if 3Z 0 2 — 27? 0 2 =0 the parabola vanishes. 
Thus one can explain the fact that the orbits a, A, b, B 
in Fig. 3 leave an open space to the right, while the 
orbits c, C, d, D leave an open space to the left. In 
case E, where 3Zo 2 = 22?o 2 (= |V) it is seen that no open 
space is left either to the right or to the left. In this case 
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Fig. 3. Orbits in the case'P = 4Q. 


776 


GEORGE CONTOPOULOS 


the angle of the velocity vector with the £ axis is 
tan- 1 (2/3) = 33°41 , 24". 

If 6 = 0 all orbits are periodic, of the form of a figure 
eight. If 6 is small (6=0.02, cases a, b, c, d) the moving 
point forms approximately figure eights, which, how- 
ever, change gradually. If Ro/Zo<0 the loops become 
more elongated, until they become very thin ; then they 
become flatter again, and so on. If R 0 /Z 0 >0 the loops 
become first more and more flat; then they begin to 
expand, and so on. In case a the orbit needs 4.4 X 10 9 yr 
(22 loops) to go from the most elongated figure eight to 
the most flat one or vice versa. We may call this time 
“the long-range half -period,” to distinguish it from the 
mean period of one loop which is given in Table V. The 
long-range half-period in case b is approximately 6X 10 9 
yr (30 loops) and in cases c and d 6.3 X 10 9 yr (32 loops). 

If b is larger (6 = 0.2, cases A, B, C, D, E) the orbits 
are again composed of loops, but these are not similar 
to closed figure eights. In this case we can only loosely 
define a long-range half-period of the order of 10-70 
X10 7 yr, comprising from 0.5 to 3.5 loops. Case E has 
the greatest number of loops in one long-range half- 
period, between 3 to 3.5. 

The characteristics of all the curves are given in 
Table IV, which is of the same form as Table I. 

The boundaries of the calculated orbits are composed 
of four arcs, except in case E. Every two consecutive 
arcs intersect at a right angle on the curve of zero 
velocity. 

The inner arc is approximately a parabola. A com- 
parison of Tables IV and V shows that for 6 = 0.02 the 
empirically found inner boundary coincides rather well 
with the parabola (75). But if 6=0.2 the coincidence is 
not good enough. Evidently higher-order terms in 6 are 
in this case important. In cases A and B the opening 
to the right is larger than indicated by the parabola 
(75), while in cases C and D the opening to the left is 
smaller than the corresponding parabola. In case E the 
parabola does not become a double line z 2 =0, but 
vanishes completely. 

Orbits with the same absolute values of the com- 
ponents of the initial velocity have the same boundaries. 


In fact, the equation of the boundary (67) depends only 
on the squares of R 0 and Z 0 and not on their sign. 
Therefore, although initially two orbits with velocities 
symmetric with respect to the z axis behave differently, 
finally the two orbits fill the same space and look very 
similar to each other. This is quite clear in cases c and d, 
and somewhat less clear in cases C, D, because the 
angular points of the boundary cannot be defined quite 
accurately in these cases. 

The maximum value of Zo 2 (3Zo 2 — 2i? 0 2 ) for a given 
value of the total velocity v 0 =(Ro 2 +Zo 2 )\ is realized 
for Z 0 2 =i'o 2 , 7£o=0; then 

z 2 =(156z> 0 2 /16<2 2 )£ 

and the opening to the right is maximum (cases a, A). 

The minimum value of Zo 2 (3Zo 2 — 2R 0 2 ) is realized 
for Zo 2 =v 0 2 /5, i? 0 2 =4z’o 2 /5. Then |Z 0 /7?o|=i, i.e., the 
angle between the velocity vector and the £ axis is 
±26 0 33'54". Then the parabola is 

z 2 =-(bv 0 2 /16Q 2 )Z 

and the opening to the left is maximum (cases c, C, d, D). 

The opening of the parabola is roughly proportional 
to 6, for the same values of \R 0 \ and | Z 0 1 . This explains 
the difference between pairs of orbits (a, A), (b,B), etc. 
Table V gives the coefficients of £ in formula (75) as 
well as the value of z given by this formula for 
£=±z>o/2()*= ±0.196. 

(b) The zero-order solution (69) gives two branches 


W-Qz 2 ) 
2v 0 Q* ’ 


(vJ—Qz-) 
2v 0 Q* 


(76) 


symmetric with respect to the £ axis. These branches 
meet at the points z=±z> 0 /(? i ; they intersect at right 
angles, because 


d£ i —Q*z d£ 2 Q*z 

dz v 0 dz vq 


(77) 


Table IV. Characteristics of the orbits represented in Fig. 3. 



£o 

zo 

Ro 

Zo 

2 F 

fmin 
(2 = 0) 

£max 
(2 = 0) 

£i 

Z 1 

£2 

z 2 

a 

0 . 

0 . 

0 . 

0.12370 

0.015301690 

-0.192 

0 . 

+0.008 

±0.391 

+0.192 

±0.077 

A 

0 . 

0 . 

0 . 

0.12370 

0.015301690 

-0.118 

0 . 

+0.110 

±0.432 

+0.177 

±0.309 

b 

0 . 

0 . 

0.05532 

0.11064 

0.015301511 

-0.194 

0 . 

+0.005 

±0.391 

+0.194 

±0.055 

B 

0 . 

0 . 

0.05532 

0.11064 

0.015301511 

-0.168 

+0.043 

+0.056 

±0.425 

+0.189 

±0.202 

c 

0 . 

0 . 

0.11064 

0.05532 

0.015301511 

0 . 

+0.195 

-0.195 

±0.020 

0 . 

±0.391 

C 

0 . 

0 . 

0.11064 

0.05532 

0.015301511 


+0.193 

-0.191 

±0.057 

—0. 

±0.391 

d 

0 . 

0 . 

-0.11064 

0.05532 

0.015301511 

0 . 

+0.195 

-0.195 

±0.020 

0 . 

±0.391 

D 

0 . 

0 . 

-0.11064 

0.05532 

0.015301511 

Z^O . 

+0.193 

-0.193 

±0.043 

+0.007 

±0.396 

E 

0 . 

0 . 

0.09581 

0.07823 

0.015299489 

-0.195 

+0.195 


£=+0.010 

z = ±0.398 
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and for z= zfci'o/<2 ! , 


For z= 0, we have 



(78) 


In first-order approximation we set 

|«>= £.+&£,• (79) 

in Eq. (67) and find after a number of operations: 

1 


>=.1 — *2 = " 


1280V 


X -840Vz 4 +0[36Eo 4 +62/? o 2 Z o 2 +51Zo 4 ]z 2 


+5r 0 2 Zo 2 (3Z 0 2 -27?„ 2 )- 


100@®»o 2 z e 


»o 2 +0z 2 

In the special case when 3Zo 2 — 27?o 2 =0, we have 


i- 


2lQz 2 +9v 0 2 - 


25Qz* 


32Qv<j t L Vq 2 -\-Qz 

For z=0, it is then E,-=0, i.e., 
Zi=±v 0 /2Qi. 

The curve of zero velocity 

4Q?+Qz 2 -2bfr 2 =v 0 2 


} 


(80) 

(81) 

(82) 

(83) 


intersects the £ axis at the points £=±0.196 and the 
z axis at the points z= ±0.391. The maximum value of 2 
is Zmax=0.391 for £=0.008 if 5=0.02, and 2 ma x=0.434 
for £=0.094 if 5=0.2. This curve is represented by a 
dashed line in Fig. 3. 

In case E, (3Z 0 2 —2Ro 2 =0), the curve of zero velocity 
and the boundary of the orbit are tangent at the points 
£i=±flo/2(?*, 2=0. 

In the general case 


£(0 = ± 


W-Qz 2 ) 

2v 0 Q* 


(84) 


v 0 55Z 0 2 (3Z 0 2 — 2R 0 2 ) 

£ (i) = ± — + . (85) 

20* 1280V 

Table V gives the values of 5H t - and £i+5Ei, £2+622 
for z=0, on the side where the boundary actually inter- 
sects the £ axis. These values coincide very well with the 
values £min(2=0) or £ ma x(z=0) in Table IV for 5=0.02, 
but less satisfactorily for 5=0.2. 

For z=dzv 0 /Q i , we have £ t =0 and 


£= b(R 0 2 - 4Z o 2 ) 2 /640 V. (85a) 


The last value is always positive (or zero) if 5> 0; this 
indicates that the point where the two branches of the 
envelope meet is to the right of the axis z. 

This point lies on the curve of zero velocity (83) to 
an accuracy of first order in 5. The accurate boundaries 
of the calculated orbits also show that the points of 
intersection of the different branches of the boundaries 
lie on the curve of zero velocity within an accuracy of 
±0.001. The value (85a) however, is not accurate, 
because if £ t =0 the approximation formula (79) is no 
more convenient. 

The points of intersection of the curves (79) and (75) 
are found by setting the value (75) in £ t - and z 2 =0 in 
Si, and solving for £. Then 


5bZ 0 2 (3Z 0 2 -2Ro 2 ) 

£= (=t ) — 1 , 

20* 1280V 


( 86 ) 


where we have the plus sign if 3Zo 2 —2R 0 2 >0 and the 
minus sign if 3Z 0 2 —2R 0 2 <0. 

If we set the value (86) in Eq. (83) we find 


55Zo 2 (3Z 0 2 -2V) 

z 2 = (±) 

32Q 5,2 v 0 


(87) 


This value is the same as (75) for £= (±)(?>o/20*). 
Therefore the intersection lies on the curve of zero 
velocity to this approximation. 


Table V. Data for the orbits represented in Fig. 3. 



3Z 0 2 — 2i? 0 2 

Coeff. of 
1 in (75) 

z for 

£ = ±0.196 

So 
11 11 

fi+KSi 

( z= 0) 

|2+iS2 

(*=0) 

Mean period 
of a loop 

a 

+0.0459 

+0.0287 

±0.075 

+0.0036 

-0.192 


20.0X10 7 yr 

A 

+0.0459 

+0.287 

±0.237 

+0.036 

-0.160 


24.4X10 7 yr 

b 

+0.0306 

+0.0153 

±0.055 

+0.0019 

-0.195 


20.0X10 7 yr 

B 

+0.0306 

+0.153 

±0.173 

+0.019 

-0.176 


23.0X10 7 yr 

c 

-0.0153 

-0.0019 

±0.019 

-0.0002 


+0.195 

19.8X10 7 yr 

C 

-0.0153 

-0.019 

±0.061 

-0.002 


+0.193 

19.1X10 7 yr 

d 

-0.0153 

-0.0019 

±0.019 

-0.0002 


+0.195 

19.8X10 7 yr 

D 

-0.0153 

-0.019 

±0.061 

-0.002 


+0.193 

19.1X10 7 yr 

E 

0. 

0. 

0. 

0. 

0. 

0. 

20 . 2 X 10 7 yr 
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Table VI. Approximate values of the third integral. 


M2 F) 

Initial 

<po 

Max 

Min 

Initial 

Max 

Min 

<po+6<pi+6 2 ^>2 

Initial Max Min 

T 

a 

2X10-* 

0. 

+ 1 -2 

0. 

+ 1.024 

+ 1.028 

+ 1.024 

+ 1.024 

+ 1.027 

+ 1.024 

5X10 9 yr 

A 

1X10- 9 

0. 

+ 19.6 

0. 

+ 10.2 

+ 16.9 

+ 10.2 

+ 10.2 

+ 15.8 

+ 10.2 

6X10 9 yr 

b 

1X10~ 9 

0. 

+0.6 

-0.5 

+0.5122 

+0.5145 

+0.5110 

0.5122 

+0.5140 

+0.5107 

6X10 9 yr 

B 

1X10- 9 

0. 

+9.3 

-5.4 

+5.12 

+8.19 

+3.29 

+ 5.12 

+6.84 

+3.05 

6X 10 9 yr 

c 

4X10" 9 

0. 

0. 

-1.1 

-0.146 

-0.145 

-0.153 

-0.146 

-0.144 

-0.149 

10X 10 9 yr 

C 

4X10 -9 

0. 

0. 

-10.0 

-1.46 

+ 1.02 

-3.13 

-1.46 

-0.07 

-3.40 

10X10 9 yr 

d 

2X10" 9 

0. 

0. 

-1.1 

-0.146 

-0.145 

-0.153 

-0.146 

-0.144 

-0.149 

1 1 X 10 9 yr 

D 

1X10" 9 

0. 

+0.7 

-9.3 

-1.46 

+ 1.38 

-2.55 

-1.46 

-0.01 

-2.81 

5X10 9 yr 

E 

2X10" 9 

0. 

+0.6 

-8.8 

-0.73 

+ 1.87 

-2.23 

-0.73 

+0.74 

-1.92 

6X10 7 yr 


At the intersection the tangent of the curve (75) is 
given by 

d$ 32QVz 

-= ( 88 ) 

dz 5bZ 0 -(3Zo 2 —2Ro 2 ) 

and the tangent of the curve (79) by 

dH Qk 

—=(±)— +*[•••]• (89) 

dz Vo 

Thus 

/d£\/da\ 32Q 5/2 » 0 z 2 

(-)( — )=(±) +••• (90) 

\dz/\ dz ) 56Zo 2 (3Z 0 2 - 2 R 0 2 ) 

and if we set the value (87) of z 2 we get 



to a zero-order approximation; we cannot go to a higher 
approximation because we have not calculated the 
higher-order terms in the expansion of z 2 . 

(c) The zero-order solution (70) can be used in the 
same way as solution (68) to give a corresponding first- 
order solution, if we set 

{=»*• 

Then Eq. (67) gives 

W*“]+2JV[-“]+*-*=0, (92) 

But now, in order to find \jz we need the terms with 6 2 
in the equation of the boundary, and these have not 
been calculated in this paper; therefore we leave this 
point open. It is probable, however, that some charac- 
teristics of the orbits, namely the fact that the boundary 
does not go through 0 in case B and similar cases, could 
be explained by the existence of another branch of the 
boundary, that corresponds to"the zero-order solution 

£(o) = 0. 

We finish with a discussion of the accuracy of the 
third integral in the calculated cases. Table VI gives 
first the accuracy with which the energy integral was 


conserved for the time interval indicated in the last 
column. Part of the inaccuracy is due to the fact that 
the calculation was eventually stopped and resumed 
later while only eight significant figures of the coordi- 
nates and velocities were available. The actual accuracy 
of the double precision calculations is about 1X10 -9 
kpc 2 (10 7 yr)- 2 in6X10 9 yr. 

The third integral is given by Eq. (20), where <po+b(fi 
is given by Eq. (48) and <p 2 is given by Eq. (54). 

The values of <?o> <?o+£<pi, and <po+6<?i+6V 2 were 
printed every 10 8 yr. Table VI gives the initial values, 
as well as the maximum and minimum values of the 
above quantities in units of 0.0001 kpc 3 (10 7 yr) -2 . 

It is seen that whenever 6=0.02 (cases a, b, c, d) the 
third integral is rather well conserved in second-order 
approximation, and even in first-order approximation, 
for time intervals of the order of 5-10X10 9 yr at least. 

The changes of the third integral to this approxima- 
tion seem to be almost periodic, therefore it is not 
probable that much larger changes will appear later. 

If 6 = 0.2 the third integral is not conserved very 
well, although its values are of the same order of magni- 
tude. This is probably because higher-order terms are 
important in this case. 

A great similarity between the maximum and mini- 
mum values of the third integral in cases c and d and a 
smaller similarity between the corresponding values in 
cases C and D is to be remarked. This verifies the fact 
that orbits with initial velocities symmetric with respect 
to the axes £ and z behave similarly in the long run. 

Case P=Q and the cases of small divisors, when 
m 2 P—n 2 Q is near zero, but not equal to zero, will be 
considered in Paper II. 
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